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Recent experimental investigations and physical modeling studies have indicated that turbulence behaviors 
within a liquid jet have considerable effects on the atomization process. This study aims to model the 
turbulence effect in the atomization process of a cylindrical liquid jet Two widely used models, the Kelvin* 
Helmholtz (KH) instability of Reitz (blob model) and the Taylor* Analogy-Breakup (TAB) secondary droplet 
breakup by O’Rourke et aL , are further extended to include turbulence effects. In the primary breakup 
model, the level of die turbulence effect on the liquid breakup depends on the characteristic scales and the 
initial flow conditions. For the secondary breakup, an additional turbulence force acted on parent drops is 
modeled and integrated into the TAB governing equation. The drop size formed from this breakup regime is 
estimated based on the energy balance before and after the breakup occurrence. This paper describes 
theoretical development of the current models, called u T-blob n and U T~TAB ”, for primary and secondary 
breakup respectively. Several assessment studies are also presented in this paper. 


Nomenclature 

a Radius of blob or parent drop 

Ai Constant (0. 1 88) 

B 0 Constant (0.61) 

Bi Constant (10.0) 

C, Constant ( B 0 /3.726B, ) 

C b Constant (1/2) 

C d Discharge coefficient of injection nozzle; constant used in equation (18) 

C D Drag coefficient of a deformed droplet used in equation (27) 

C F Constant (1/3) 

Cjd Constant (10) 

C k Constant (8) 

C t Empirical constant involving turbulence force 

Ct Weighting parameter associated with turbulence motion 

Cw Weighting parameter associated with surface motion 

c M Turbulence constant (0.09) 

C E Turbulence constant ( 1 .92) 

Maximum possible diameter of product drop 
d^ Minmum possible diameter of product drop 
D Diameter of injection nozzle 

E Term associated with energy 

F Aerodynamic force acting on the drop 

F t Force associated with turbulence motion in a liquid droplet 

k Turbulence kinetic energy per unit mas in the k — 8 model 
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kc Coefficient parameter ( o.75 / L t ) 

k( Tuibulent kinetic energy 

kta Kinetic energy associated with turbulence motion 

kfcw Kinetic energy associated with surface wave motion 

K Constant (10/3) 

K* Loss coefficient due to nozzle inlet geometry 

K e Turbulence constant (0.27) 

L Length of injection nozzle 

Ljet Liquid jet intact length 

Lw Characteristic length scale associated with surface motion 

Lt Characteristic length scale associated with turbulence 

m Mass of a liquid drop 

nibr Stripped mass of "blob 7 ’ drop during primary atomization process 

r p Radius of a parent drop used in the secondary breakup formulations 

Radius of a product drop used in the primary breakup formulations 
r t Radial length scale associated with turbulence motion 

r w Radial length scale associated with surface wave motion 

Re Reynolds number ( 2p g Wi^ J\i % ) 

s Contraction area ratio of injection nozzle 

SMD Sauter-mean diameter of product drop 

SMR Sauter-mean radius of product drop 

t Time 

t* Characteristic breakup time ( 2r p ^p,/p g /W o ) 

x Axial coordinate on the liquid injection direction 

U Liquid jet injection velocity 

v 1 st normal velocity component of product drop with respect to the path of parent drop 

V blob or parent drop velocity 

w 2 nd normal velocity component of product drop referenced to the path of parent drop 

W Relative drop velocity with respect to local gas velocity 

y Non-dimensional parameter of distortion displacement (£/C b r p ) 

Greek Letters 

AE Change in energy of parent drop 
A Wave length of fastest growing wave 

p Density 

G Surface tension coefficient 

£ Dissipation rate of turbulence kinetic energy per unit mass 

£ t Dissipation rate of turbulence kinetic energy 

X Time associated with surface wave motion 

T w Characteristic time scale associated with surface wave motion ( a/ ) 

T t Characteristic time scale associated with turbulence 

0 Jet divergence or spray angle 

p Molecular viscosity 

CO Angular frequency 

Q Radial cross-section change of drop from its equilibrium position 

Q Maximum wave growth rate 

Subscripts 

af Parameter after mass stripping from blob drop in primary breakup 

be Parameter before mass stripping from blob drop in primary breakup 

* Parameter defined at time of breakup 
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Initial value 
Value of ambient gas 

Parameter associated with motion kinetic energy 
Value of liquid phase 

Parameter associated with parent drop 
Parameter associated with product drop 


t Parameter associated with turbulence 

to Parameter associated with turbulent kinetic energy 

^ Parameter associated with surface tension energy 

w Parameter associated with surface wave motion 

Superscripts 
o Initial value 

. First derivative with respect to time 

.. Secondary derivative with respect to time 


I. INTRODUCTION 


T HE transformation of a liquid body into droplet sprays in a gaseous surrounding is of great importance in many 
industrial processes. In the liquid fuel combustion systems, such as rocket engines, diesel engines, gas turbines, 
and industrial furnaces, the combustion efficiency and chemical reaction behavior are primarily dependent on the 
effectiveness of the liquid body broken up into sprays. Smaller drop sizes generated from the spray devices increase 
the specific surface area of the fuel and thereby achieve high rates of mixing and evaporation. On the other hand, 
when fuel is mixed and then reacts rapidly near the injector, the injector exit surface could be overheated 
inadvertently. These phenomena have been observed in many liquid rocket engines during their hardware 
development phase. Therefore, understanding and adequately predicting this physical breakup process leads to 
designing better spray devices for these various applications. 

The mechanisms of atomizing the liquid jet are complex. The jet inertial and aerodynamic forces along with the 
surface tension certainly play a role in the breakup process. In addition, the turbulence behavior inside the liquid jet 
also contributes to the jet disintegration. Often the geometrical sharpness of the injection nozzle inlet, along with 
appropriate flow conditions, can create cavitations inside the nozzle. The collapse of this cavitation can generate a 
flow fluctuation, leading to a more aggressive disintegration of the liquid jet. Previous studies [1-4] showed that the 
breakup length of the liquid core was significantly shortened when the liquid jet changed from laminar to turbulent 
flow. More recent studies of relatively large-size liquid jets [5-8] provided droplet data such as size, velocity, 
fluctuating quantities and reported a more complete mechanistic approach for the primary atomization regime, and 
indicated that the turbulence developed inside the jet column, starting at the nozzle exit, remained dominant and 
became a main contributor in the spray development. The authors believed that this turbulence characteristic played 
a primaiy role on the liquid stripping near the injector exit Based on the data collected from their experiments, the 
authors were able to establish correlations between the turbulence fluctuation quantities and the breakup drop size 
and the breakup length of the liquid jet. 

In the modeling arena, as to the convenience of implementing into computational fluid dynamic (CFD) methods 
and the wide utilization by analysts, the two most noticeable atomization models are the KH instability model of 
Reitz [9] and the TAB of O’Rourke et al. [10]. Reitz derived the KH model, also known as the “blob” model or 
‘‘stripping rate” model, which described the primaiy breakup entirely in terms of the surface wave perturbation, in 
which the viscosity, surface tension, and aerodynamic forces were the contributing factors. On the other hand, the 
TAB model was based on an analogy between an oscillating, distorting droplet and a spring-mass system. In this 
model the drop distortion was driven by the force interaction among the externa! aerodynamic, surface tension, and 
viscous damping of the droplet liquid. The TAB model is suitable for predicting the secondary breakup regime. 
Recently, Tanner [11] examined data from liquid jet breakup experiments at high pressure conditions and found that 
an intact liquid core is either broken into various liquid shapes or drop sizes shortly after the injection exit port. 
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Hence, the author extended the TAB model and developed a so called enhanced TAB (ETAB) model, to include the 
primary breakup regime. In his recent Cascade Atomization Breakup (CAB) model for high-pressure liquid jets. 
Tanner [12] also accounted for the droplet-surface stripping near the injector exit in his ETAB model with a power 
law drop-size distribution. Though these models provide reasonable predictions of the liquid atomization in some 
aspects, they do not account for the liquid turbulence motion observed in certain sprays as previously described. In 
a study of the primary diesel fuel atomization, Nishimura et al [13] developed a phenomenological cavitation model 
in which the primary breakup is governed by the turbulence mechanism. The turbulence energy is formed by the 
bubble collapse and the fluid turbulence motion. To consider the nozzle exit turbulence conditions in the modeling 
of diesel sprays. Huh et al [14] proposed a scheme taking into account two independent mechanisms, wave growth 
and turbulence in the atomization process. The turbulence is characterized partially by the injection nozzle 
geometry while the wave growth is derived hom the KH instability theory. With this attempt a logical framework 
of coupling the flow inside the injection nozzle to the jet breakup can be achieved. The rationale for this approach is 
based on an order-of-magnitude analysis of the flow dynamic breakup mechanisms and the ones associated with the 
turbulence. This analysis concludes that the gas inertia force and turbulent stress are the two main forces in the 
atomization of typical diesel engine injectors. In this model, the breakup rate for the parent drops is set to be 
proportional to the length scale to time scale ratio. It is hypothesized that the turbulence length scale is primarily 
dominant in the primary breakup while the wave length scale is important for the secondary droplet breakup. On the 
other hand, the time scale is a linear sum of the turbulence and wave growth time scales. 

The aforementioned models were formulated with semi-empiricism. The associated coefficients had to be 
determined with reference to experimental data. Another class of atomization models was derived based on the 
conservation laws of mass, momentum, and energy. Several efforts of using sophisticated numerical approaches 
have been performed in modeling the detailed turbulent flow fields in the liquid and gas during the atomization 
process. Klein et al [15] and De Viiliers et al [16] computed the two-phase flows using the large eddy simulation 
(LES) method. For the first time, De Viiliers et al. were able to apply such a numerical technique to resolve the jet 
atomization under more realistic operating conditions found in Diesel engines. Leboissetier et al [17] also 
attempted to simulate a liquid spray with a multi-dimensional pseudo-direct numerical simulation (DNS) method. 
These methods have a real potential of providing a complete physical description of the liquid jet breakup with 
minimum assumptions; however, they require small computational time steps and fine grids across the entire jet 
domain for their simulations. Especially the LES and DNS techniques may need submicron spatial elements in size 
and pico-second in time steps, to properly predict the atomizing sprays at the high-velocity injection conditions. 
Consequently, grid mesh size and die considered physical domain must be taken care of so that the computational 
time and memory storage requirements can be manageable. At the present time they are still too expensive and 
generally impractical in terms of computational time and power requirements for engineering calculation 
applications. Hence, the engineering analysis and design of the liquid spray devices still must rely on 
phenomenological engineering models. 

As indicated earlier, the results of experiments suggest that the turbulence motion within the liquid may alter the 
intact jet core as well as the production of drops, which have a different size in comparison with the droplet breakup 
without the turbulence consideration. Thus, an accurate modeling of such turbulence effects on the atomization 
process would significantly enhance the prediction capability of existing physical models. The purpose of this 
research aims to model the turbulence effects on the primary and secondary atomization processes. In the course of 
this study, terms accounting for the turbulence motion within a liquid are developed, based on a phenomenological 
point of view. These terms are appropriately supplemented to the two classical atomization models: HK primary 
breakup [9] and TAB secondary breakup [10], respectively. In the primary atomization model, the level of the 
turbulence effects on the mass stripping of die blob drops and on the product drop size is represented by the 
characteristic time and length scales and the initial kinetic energy. This treatment offers contributions of individual 
physical phenomena on the liquid breakup. For the secondary breakup an additional turbulence effect, acted on the 
parent drops, is modeled and integrated into the TAB governing equation. This turbulence term is referenced to the 
dissipation rate of the turbulent kinetic energy and the deformation rate of the parent drop distortion displacement. 
The results of the new models and the existing ones are compared and appraised using available experimental data. 
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EL MODEL DEVELOPMENT 


A. T-blob Model 

The original formulation of the primary breakup blob model is based on the stability analysis performed by Reitz 
[9]. The main contribution of this research effort is to incorporate the turbulence effect in modeling of the liquid jet 
breakup. Terms and parameters associated with the turbulence behavior are derived and implemented appropriately 
to the existing blob model. Hence, it is appropriate to denote the present model as "T-blob” model. In order to 
account for the both effects of the surface wave perturbation and the turbulence motion, the resulting formulation 
would include the combination of these two phenomena in die breakup process. The length scale and time scale 
associated with the primary breakup are comprised of the ones described in the blob model that represents the 
surface wave instability and of the turbulence behavior following die approach of Huh et al. [14]. To derive the new 
model, it is proposed that the rate of change in the parent drop radius has an extra term associated with turbulence 
effect in the form of 
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n 
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,when 


r p <^a. 


( 1 ) 


where 


L W =A 
x w =a/AO 
x = 3.726B,a/A£2 

C, = — B ° 

* 3.726B, 

The characteristic length scale L* and time scale T w are associated with the droplet surface wave instability. These 
two scales along T and C a are formulated from the u blob ” model. The turbulence characteristic length scale L, and 
time scales T t can be derived using the analytical solution of the k — £ turbulence model. 

T t = x t ° + 0.0828t (2) 


L,=L“ 


1 + 


0.0828t 


x 0.457 


( 3 ) 


The time t is counted from the time at which the parent drop leaves the injection nozzle exit. The initial turbulence 
length scale L° t and time scale are evaluated from the initial turbulent kinetic energy k° t and its corresponding 

dissipation rate E° t at the injector exit 

L"-cliC 

0 k? 

X « =C hT5-’ 


where 


<^=0.09. 


The initial turbulent kinetic energy k° and its corresponding dissipation rate can be approximated as 
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( 5 ) 


e° = K — 

' 2L 




with K £ = 0.27 - The velocity U is the liquid velocity at the injection nozzle, which has the length, L, and the 

diameter, D. The discharge coefficient, the loss coefficient due to the nozzle entrance sharpness, and the 
downstream-to-upstream contraction area ratio of the injection nozzle are represented by C d , K*, and s, respectively. 
Detailed derivations of these turbulence scales can be found in [14, 21]. 

The third term on the right hand side of equation (1) represents the contribution of the turbulence effect on 
the stripping rate of the parent drop. This turbulence term would accelerate the reduction in the parent drop radius 
with time and tire modeled effect is consistent with tbe experimental observations [3, 18-20], in which the liquid jet 
has a short intact core length when the turbulence appears in the liquid jet due to the flow at a high Reynolds number 
or cavitation. It should be noted that for the non-tuibulence case, equation (I) would become the original blob 
model of Reitz, which only retains Che effect of the surface wave motion on the primary breakup. Similar to the 
criterion set by Reitz, die parent drop would no longer strip its mass to create the product drop when the parent drop 
radius a is less than the radius r p of the product droplet 

Studies of turbulent liquid jets [5, 6] indicated that the drops produced from the onset breakup regime have their 
sizes in the same order-of-magnitude of the turbulence length scale. On the other hand, Reitz set the radius of the 
product drop to be proportional to the wave length associated with the fastest growth rate. The above relations 
suggest that it is reasonable to formulate the radius of the product drops with the characteristic length scales for both 
the wave perturbation and turbulence phenomenon. For the present model, the reciprocal of the product drop radius 
is expressed by the sum of the reciprocals of the length scales associated with the surface wave instability and the 
turbulence motion with the inclusion of the respective weighting factors as follows 


2 

r P 



r w r , 


(l-e,)r,+c t r w 


( 6 ) 


Again, the radius of the product drop is r p . The radius r w associated with the wave motion can be determined 

from the Reitz model. In die present phenomenological description, the kinetic energies associated to the surface 
wave and turbulence characteristics are used to weight both effects. The consideration of the kinetic energy level for 
this weighting is based on a reason of which the phenomenon of the larger kinetic energy motion would have a 

stronger influence in the liquid jet breakup process. Hence, the weighting coefficients C w and C t are determined by 


the kinetic energy ratio of the turbulence motion and wave perturbation. The value of r t is estimated from a 

probability density function (PDF), which was proportional to the ratio of the turbulence energy spectrum and the 
atomization time scale [14]. The notion for this representation was that an eddy motion with larger turbulent kinetic 
energy and a shorter atomization time most likely caused the drop breakup to occur more frequently than other ones 
containing the lower energy level and the longer atomization time. By assigning the wave number of the turbulence 
energy spectrum as the inverse of the product drop diameter, it can be shown in [21] that 
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where k e = 0.75/L t , = maximum possible diameter of the product drop, and 6^= smallest possible diameter 

of the product drop. The weighting coefficient C| is formulated with the kinetic energy terms related to the surface 
wave and flow turbulence levels as 

c, = * , (8) 


dx 


Eb + E kw 


where 
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Since the weighting coefficients represent degrees of the contributions of a particular physical phenomenon to 
the overall one, it is logical to assign a unity relation between them as C w + C t = 1 . With this relation, it is easy to 
recognize that for the constraint of E ta » E^ , r p would become equal to r t , and vice versa r p would approach r w if 
E fa <3C E,^ . This ensures that the condition with a larger kinetic energy would have a more significant effect on the 
atomization process. 

Along with the inclusion of the turbulence effect on the primary atomization process, the subject phenomenon 
also is considered in the secondary droplet breakup model, which is discussed in the next section. At any rate this 
new model requires the initial velocity fluctuation quantity of the product drops right after their formation. This 
quantity can be obtained by examining the energy conservation during the primary breakup process. Denoting 
and as the radii of the parent drop before and after its liquid stripping, respectively, the amount of the mass 
mbr stripped from the parent drop is formulated as 

4 


“br=7^P.(a 3 be-a!f)- 


(9) 


During this breakup process, die change in die individual energy forms of the parent drop is estimated as follows: 


Energy due to surface wave motion: 
Surface tension energy: 

Turbulent kinetic energy: 

Kinetic Energy of motion: 


(AE. )„=«*( Afi) 

( AE -f)_ =2,ro ( , £- , 2) 


V T t 
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/par U1 ^ 

For the product drop, the energy associated with the surface distortion is negligible as compared to the other energy 
forms, since the initial disturbance is small. Hence, the energy of the product drops is composed of the surface 
tension, turbulence, and kinetic motion as described below. 


Surface tension energy 


Kinetic Energy of motion 
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where and are the velocity and number of product drops, respectively. They are evaluated as 

( 10 ) 


= Vv 2 +w 2 


prod 

N - faV 
^ 4*prf 


( 11 ) 


For the velocity, the product drops would cany the same velocity V of the parent along with two additional 
normal velocity components, v and w, in reference to the trajectory of the parent drop. They are represented by 

v = | V| tan (0/2) sin <j) 
w =|v|tan(0/2)cos(j>, 

where 
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tan(e/2) = A,Afi/U. 

A random parameter $ is chosen on the interval (o,2u) - According to Reitz et al. [40], the constant A x depends on 

nozzle design. Its value is set equal to A|= 0. 1 88 for the present study. 

By equating the change in die energies of the parent drop with the energy of the product drops, the turbulent 
kinetic energy for the product drop is written as 


< E -),»d = 


(AE.) +(AE„, ) +(AE„ ) +(AE t ) 


/par 


N, 


prod 
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ae_U+(ae,) i 


prod 
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( 12 ) 


This quantity is used to estimate die initial velocity fluctuation for the secondary breakup process, which is 
discussed next 


B. T-TAB Model 

Studies [7, 8] of the turbulence effect on the liquid breakup suggest that the turbulence motion tends to weaken 
the surface tension force. In fact, it is well recognized that this surface tension keeps the liquid drop from being tom 
off, while the turbulence within the drop would promote the droplet disintegration process. To account for this 
behavior, a term F t , representing a force associated with the effect of the turbulence on the droplet breakup, is 
introduced to the original TAB equation 

m£ = F+F t -k (13) 


It is assumed that the turbulence force F t is formed by a portion of the internal liquid turbulence energy. Another 
portion of subject energy decays through a dissipation process. This proposition is based on the experience that the 
turbulence motion behaves like a force participating in the droplet surface deformation and accelerating its surface 
displacement Hence, it is proposed that die product of the turbulence force and the deformation rate of the droplet 
surface displacement must be related to the dissipation rate of the turbulence energy during this breakup process. 
This relationship is formulated as 

F,4 = C.me (14) 

The dissipation rate of the turbulent kinetic energy per unit mass of individual drops is £ , and an empirical constant 
representing the proportionality of subject relationship is C t . 

It should be noted that without additional turbulence generation the internal turbulent kinetic energy would decay 
with time as the drops travel downstream. With an assumption of the turbulence within the drop being homogenous 
and isotropic, using the k- £ model [22], k can be expressed as a function of time: 
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(16) 


The turbulence force, shown in equation (14), can be rearranged and written as 

- L =C,£ 0 
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( 17 ) 
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When the distortion displacement t, is non-dimensionlized by y = C b r p , equation ( 1 3) becomes 


im £L£iK+ c 


' Ctf L 


y'-C k 


-^ry-Ca-^y 

p. r p p. r p 


( 18 ) 


<- b p. r p 

Except for the second term on the right hand side, equation (18) is the same as the governing aquation of the TAB 
model. To estimate the post-breakup drop size, a similar method as used in the TAB model is employed for the 
calculation of Sauter mean radius. 


SMR = 


4ot p 3 o 


r? ft 5 - 2 

E p> r -^ r P Piy 


The total energy of the original drop prior to the breakup is , and it is determined from equation the original 
TAB formulation 


E par = E surf + E osc = K-ftt 3 ' + ©V) - 


(19) 


The secondary breakup process is considered for all drops created from either the primary breakup or previous 
secondary droplet breakup and requires initial values of the turbulent kinetic energy and its dissipation rate. When a 
drop is created from the secondary atomization, it is assumed that the turbulence quantities of the parent drop are 

preserved and distributed evenly to its children drops. The turbulent kinetic energy k bt and corresponding 

dissipation rate 8^ of the parent drop at the breakup time are determined from equations (15) and (16). The time and 

initial turbulence values used in these equations are referenced to the time when the considered drop is formed. The 
formulations employed for this calculation are 




(k.f 


0 C e)tl* + 0 C o) 
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P = H. 

bt lf C ‘ 
K o 
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(k r 


(l-c e )t bt +(k 0 )' 
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( 20 ) 


( 21 ) 


Again, t&, ko, and £ 0 are the breakup time and the initial turbulence quantities of the parent drop. Based on the 

conservation of mass, the initial turbulent kinetic energy and dissipation rate of the new drop are simply formulated 
as 


k =— L k 

0 SMR 3 


and 8 — — 

° SMR 


3 £ bt * 


( 22 ) 


It should be emphasized that ko ande o are denoted for the initial turbulence quantities of a corresponding drop type. 

The parameters in equations (20) and (21) belong to the parent drops, while they are for the product drops when 
used in equation (22). 

When a drop of interest is directly produced from the primary atomization, its initial turbulence energy is 
determined from the energy conservation of the blob drop and its product drops at the breakup time. Hie 

formulation of the initial turbulent kinetic energy for the product drop, denoted as(E tur , has been derived and 

presented already in equation (12). Hence, the initial turbulent kinetic energy k« per unit mass for the drop is 
represented by 

k.=-4-(E»*)„ 


4mt p Pi 


/prod * 


The initial dissipation rate of the turbulence energy £ 0 is then determined from the relationship 
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With the assumption that L t = 2r p , £ 0 can be written as 


e o C M 


2r„ 


It should be emphasized that these two relationships are used for estimating the turbulence values of the new drops 
for the primary breakup only. 

In an atomization process without the consideration of the droplet coalescence, the liquid drop breaks up into 
smaller drops. As a drop travels downstream its size reduces to the smallest possible dimension, and the turbulence 
activity within the drop also decays. Subsequently, it would lead to conditions in which the turbulence approaches 
its final energy cascade regime. Here, the molecular viscosity is effective in dissipating the turbulent kinetic energy, 
and the turbulence scales also reduce to the order of the Kolmogorov scale magnitude [23], Hence, a criterion for 
eliminating the turbulence effect is postulated as follows 

3 I 1 /* 

Hi 


SMD<4 


pf £ o 


(23) 


The bracketed term with its exponent represents the Kolmogorov length scale and is based on the liquid 
properties and the initial dissipation rate of the turbulence energy of the new drop immediately after the droplet 
breakup. When the above expression is satisfied, the turbulence effect is no longer valid and the T-TAB model 
would become the original TAB model. As seen from the results of all test cases in this study, the new drop size 
reaches the Kolmogorov scale in about two or three cycles of breakup steps. 


III. RESULTS and DISCUSSION 

A. Primary Breakup Assessment (T-blob) 

Several experiments were selected to assess the current model. The flow conditions of the cases are summarized 
in Table 1. The turbulence scales, initial turbulent kinetic energy, and its dissipation rate required in T-blob model 
are estimated from equations (4) and (5) by knowing die injection nozzle configuration and its flow conditions. The 
values of die kinetic energy and the dissipation rate are listed in Table 1 . It should be noted that for some test cases 
the measured data to calculate the initial turbulence values are not entirely reported. Subsequently, typical values 
are estimated and utilized in the calculation. Details of determining these initial values can be found in [21]. 

We first investigate the relative effects by comparing the mass stripping rate due to turbulence and surface wave 
as described in equation (1). The values of these terms for test cases H-l and H-3 (see Table 1) are shown in Figure 
1. We first investigate the relative effects by comparing the mass stripping rate due to turbulence and surface wave 
as described in equation (1). The values of these terms for test cases H-l and H-3 (see Table 1) are shown in Figure 
1 . In this figure the values of these terms are plotted against the relative life time of the parent drop. This parameter 
is non-dimensionalized by the injection velocity and injection nozzle diameter. The results for the considered cases 
indicate that the surface wave perturbation has a considerable effect on the reduction rate of the parent drop size. 
The value of the wave motion term in the test case H-3 is approximately two orders of magnitude greater than the 
one of the turbulence term. However, when the injection velocity increases from 86.4 m/sec (case H-3) to 102.0 
m/sec (case H-l) the gap of these two values becomes smaller. This suggests that the value of the turbulence term 
rises at a faster rate than the one of the wave motion term, when the injection velocity increases. The reduction rate 
of the parent drop decreases with the increase of the velocity. Subsequently, the primary atomization process is 
elongated for the case of higher injection velocity. The reflection of the aforementioned effects can also be seen in 
Figure 2, where the parent drops predicted by the HK (blob) and T-blob models for test cases HI through Y-3, are 
plotted. 
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Table 1. Test cases and measured data used in the computations 


Case 

H-l 

H-2 

H-3 

Y-l 

Y~2 

Y-3 

K 

s 

Nozzle Diameter (mm) 


0.3 



0.213 


0.24 

0.15 

Ambient Gas 


Nitrogen 


Carbon Dioxide 

Nitrogen 

Nitrogen 

Ambient Pressure (MPa) 

U 

3.0 

5.0 

4.5 

2.5 

0.5 

2.17 

1.5 

Ambient Temperature (K) 


298 



298 


298 

289 

Density (kg/m 3 ) 

12.36 

33.70 

56.17 

72.61 

40.34 

8.07 

24.51 

16.84 

Liquid Fuel 


Diesel Fuel 



Diesel Fuel 


Diesel Fuel 

Diesel Fuel 

Density (kg/m 3 ) 


840 



840 


840 

840 

Viscosity (Kg/m' s) 


2.9x1 O' 3 



5.0x1 O’ 3 


5.0x1 O’ 3 
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Figure 1 . Values of the wave motion and turbulence terms on the reduction rate 
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It is evident from Figure 1 that the reduction rate in the parent drop size is nearly constant with time. This trend 
is due to the fact that the term, associated with the surface wave motion, is determined from the wave length and its 
fast growth rate, which are not a function of time. On the other hand, the turbulence term is time dependent. 
Combining these two equations when calculating the turbulence term, however, only offers a slight decrease in its 
value with time. In addition, the end of the parent drop breakup occurs in a relatively short time. Hence, the results 
show the rate of change in the parent drop size is almost constant. Therefore, the curves of the parent drop 
diameters as shown in Figure 2 are nearly straight downward. In all test cases, the T-blob model predicts the 
completion of the parent drop breakup slightly earlier than the prediction from the blob model. 

It should be pointed out that the T-blob model only describes the primary breakup. In the actual measurement, it 
is, however, difficult to separate other physical phenomena, such as the secondary droplet breakup, etc. Therefore, 
only die intact core length of the liquid jet can be legitimately used to compare the measurement and prediction. So, 
an available correlation [26] of the intact length, which has been widely used in literature, is used to compare with 
the present prediction. The correlation is shown below 



The intact length is denoted by and the constant has a value of 10.0. Due to the nature of the implementation 
of the blob and T-blob models, it is reasonable to define the intact length of the liquid jet as the traveling distance of 
the injected blob droplets during this process. The non-dimensional intact length L j ct /D , predicted by both 

models, plotted against the square root of the density ratio ^p,/p g , are shown in Figure 3. 
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Figure 2. Rate of change in parent drop size predicted by KH (blob) and T-blob model 
(a) Test cases H-l, H-2, H-3 (b) Test cases Y-l, Y-2, Y-3 
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Figure 3. Comparison of predicted and correlation intact lengths of the liquid jet 


The predicted results from the blob and T-blob models reasonably agree with the correlation curve. In general, 
the data points from the T-blob model are closer to the correlation line in comparison with the ones from the blob 
model. Furthermore, the T-blob model predicts a shorter intact length than the blob model does. This prediction is 
consistent with the measured data trends observed by other authors [3, 26]. 
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Figure 4. Parameters used to determine product drop size 

(a) Kinetic energies and weighting coefficient; 

(b) Radial length scales and weighting coefficient 


In the present model, both the liquid jet surface wave perturbation and the turbulence motion play a part in 
forming the product drops. Their drop size, as shown in equation (6), is composed of the two radial length scales r w 
and r t . They are weighted by the kinetic energies of the respective phenomena when used to determine the radius of 
the product drops. The values of subject parameters for the test case H-3 are plotted in Figure 4. The curve showing 
the kinetic energy of the surface wave motion in Figure 4(a) is every much constant with time, because the surface 
wave length and its corresponding growth rate are not a function of time. On the other hand, the kinetic energy 
associated with the turbulence is initially at a high level due to the turbulence developed at the injector nozzle exit. 
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However, its value rapidly decreases through the dissipation process and approaches the level of the surface wave 
motion. This trend has led to define the weighting coefficient curve (see equation (8)), which starts at the value near 
unity and then quickly reduces its value to one half at the end of the product drop generation. Figure 4 (b) presents 
the radial length scales formulated from the two considered motions. Using the same previous argument of the 
surface wave length being constant, it is easy to recognize that the radial length scale r w associated with this motion 
is more or less constant with time. In contrast, the value of r t involving the turbulence increases with time since the 
turbulence length scale Lt increases with time and the radial length r t is a function of L*. When applying the 
weighting coefficient Q (its curve is also shown in Figure 4 (b)) in the determination of the product drop size, the 
resultant value of r p initially rises and then gradually drops due to the reduction of Q. This trend can be interpreted 
such that the high turbulence intensity at the initial primary breakup stage controls the formation of the product drop 
size. As the parent drop travels downstream the turbulence dissipates. Then, this drop formation process is 
gradually dominated by the surface wave perturbation. For the test case H-3, r w determined from the Reitz model 
has a value less than 1 Jim while r t formulated by Huh et al . [14] is approximately 20 Jim. Due to the ignorance 

of the drop collision and coalescence effects, the value of r t is more representative of the drop size found in the 
measured data. 

This same value is also predicted by the T-blob model at the initial product drop formation. However, the drop 
size continues to decrease to the level predicted by the blob model. As the parent drop strips its mass and 
consequently reduces its size when traveling downstream without coalescence, it is logical to expect that the product 
drop size should decrease also in this process. It should be noted that the same trends of all parameters describing 
test case H-l are also observed in all the test cases shown in Table 1 . The prediction of the product drop sizes for 
test cases H-l to Y-3 is displayed in Figure 5. 

The results indicate that when liquid is injected at a higher velocity (test cases H-l and Y-l) product drops are 
initially formed with a larger size and this product drop formation process is taken for a longer time in comparison 
to the process at the lower injection velocity (test cases H-3 and Y-3). Similar to the observation from Figure 4, the 
results shown in Figure 5 reveal that the product drop size predicted from the T-blob model is approximately one 
order magnitude larger than the one from the blob model. Again, by ignoring other effects, the drop size estimated 
from T-blob is consistent with the experimental data. 
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Figure 5. Product drop size predicted by KH and T-blob models 
(a) Cases H-l, H-2, H-3, (b) Cases Y-l,Y-2,Y-3 


As the product drops from the primary liquid jet breakup continue traveling downstream these drops are exposed to 
external as well as internal forces. Subsequently, they may undergo additional breakup cycles. This phenomenon is 
driven by the secondary breakup mechanism. The modeling assessment of which will be discussed in the following 
section. 
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B. Secondary Droplet Breakup Model (T-TAB) 

To assess the secondary droplet breakup model, the breaking up of an isolated liquid droplet traveling in a 
gaseous medium was numerically calculated [21]. The model is computed with a temporal variation of relative 
velocity between the drop and gas medium. Hence, the force terms in the governing equation of the drop breakup 
process can be predicted in a more realistic manner. 

We first evaluate the turbulence force constant described in equation (14). The coefficient Q represents the 
proportionality between the dissipation rate of the turbulent kinetic energy and the turbulence energy participating in 
the droplet deformation. For the present study, the value of Q is estimated from measurements of drop breakup 
provided by Chou et al. [27]. In their shock tube experiment, the breakup of a single liquid drop in a flow behind a 
moving shock wave was investigated. The gas in the driven section flows at a velocity of approximately 80.8 m/s 
with a density of 1 .48 kg/m 3 . The average product drop SMR of the entire secondary breakup process and the parent 
drop axial velocity Up are collected and reported in forms of the following correlations 



The characteristic breakup time t* is defined as 2l v ^Jp ] Jp g / W G . The time t is referenced to the time at which 

the drop is introduced to the flow field, while is the time when the secondary breakup is complete. The results 
from subject experiment suggested that tt*/t* was approximately equal to 5.5. Therefore, this value is used in 
equation (25) for estimating the SMR of the product (hops. In the present study, the relative velocity of the parent 
drop is computed from equation (26) for the secondary breakup models (TAB and T-TAB). Due to the lack of 
measured droplet turbulence quantities, the initial turbulent kinetic energy ko used in the T-TAB equation is 
estimated from the parent drop velocity. Based on numerical simulations of the primary atomization for cases listed 
in Table 1 , the ratio of the fluctuation velocity within the drop and the drop velocity ranges from 0.09 to 0. 1 1 [21 ]. 
The ratio of 0.1 is selected for determining the initial droplet turbulence energy. Since Chou et al . [27] and Liang et 
al. [28] indicated that the drop starts the breakup at tit =1 .5, the initial drop velocity is predicted for this particular 
time. The initial dissipation rate £ 0 required in T-TAB equation was also estimated based on the discussion in 
Section II-B. A single water drop of 590 (im in diameter has been analyzed using the density, viscosity, and surface 
tension constants of 997 kg/m 3 , 8.94X10 -4 kg/ms, and 70.8xlCT 3 N/m, respectively. The product drop size for this 
case at various C, values is plotted in Figure 6. Since the correlation of the experimental data is not related to Q at 
all, its value shown in this figure at SMR/r p =0.135 remains constant. Note that the TAB model does not carry the 
turbulence term; therefore, the results from this model would also not vary with Q. Furthermore, TAB predicts the 
product drop size at SMR/r p =0. 145, which is roughly 7.5% larger than the measured value. The solid-line curve in 
Figure 6 displays the results of the present model with a variation of Q. As would be expected, the T-TAB model 
reproduces the same results of the TAB model when Q = 0. With an increase in the value of Q, T-TAB predicts a 
smaller product drop size. In other words, the results of the present model suggest that the secondary atomization 
process would create a smaller drop size when stronger turbulence exists within the parent drop. From the plot, it is 
evident that the predicted SMR of the product drop matches the measured data for C t =0.19. Consequently, this 
value is selected for the turbulence constant used in this study. 
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Figure 6. Variation of product drop size due to turbulence force coefficient Q 


To further assess the T-TAB model, the relative velocity between the drop and the gaseous environment is 
calculated from following relationship: 


dW_ 3C D p gw2 
dt 8r p p, 


(27) 


The drag coefficient C D is determined from a model offered by Liu et al [29]. In this model, the drag of a deformed 
droplet is composed of two parts. The first part is the drag of a rigid sphere, while the second part involves drag 
associated with the drop surface deformation. The relationship is represented as 


c D =< 


IK Rc4n )( ,+2 ' 632y ) 

0.424(1 + 2.632y) 


Re <1000 
Re >1000 


(28) 


The Reynolds number Re is defined as 2p g Wr p j |X g . The droplet distortion displacement y is calculated from 

equation (18). When the value of Co is computed from a previous time step, equation (27) becomes an ordinary 
differential equation and can be solved easily [21]. Since the initial turbulence quantities are required for the T-TAB 
model it is worthwhile to conduct a sensitivity study of the turbulent kinetic energy ko on the predicted results. 

Variations of the product drop size and drop breakup time at different k 0 /.5mW 0 2 values for test case H-3 in 


Table 1 are plotted in Figure 7. 

The x-coordinate represents the ratio of the initial turbulence energy ko and droplet motion kinetic energy, 
1 j 2mW o 2 . The curves of both product drop size and breakup time show that these quantities are decreasing almost 

linearly with an increase of the initial turbulence energy, except for k, ~ 0. These results suggest that the level of the 
initial turbulence within the liquid has a strong effect on the droplet breakup. It is also of interest to note that the 
secondary atomization process generates smaller drops with the inclusion of turbulence. Also, it takes a shorter time 
to break up the parent drop when turbulence is considered. When implementing the secondary atomization model 
with the CFD code, the initial turbulence quantities are estimated from the process of the parent drop formation. As 

discussed in [21], CFD results showed k 0 /.5mW 2 ranges from 0.09 to 0.11 for test case H-3 when the drop is 
formed during the primary atomization. Hence, it is reasonable to assign k 0 /.5mW 0 2 = 0.1 for the current study 
here. 
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In the T-TAB model, the droplet distortion motion is driven by the interaction of aerodynamic, turbulence, 
surface tension, and viscosity events. The budget of these quantities with time is shown in Figure 8. 



Figure 7. Product drop size and breakup time of various initial turbulent kinetic energy 



Figure 8. Forces generating droplet distortion based on T-TAB 

This plot depicts that the aerodynamic force has a strong effect on the droplet breakup process. The flow 
conditions for the test case H-3 are at a Reynolds number ( 2 p g r p W o /p g ) and Weber number ( 2p g r p W Q 2 ja ) of 

3.8X10 3 and 80.7, respectively. It should be noted that the time derivative of the distortion displacement y, which 
can be seen later in Figure 9, starts at a small value and increases with time. Since the turbulence force is inversely 
proportional to the y-derivative, this force consequently is very high at the initial time and rapidly changes to a 
smaller value. In contrast, the surface tension and viscosity forces are initially small and change to large values with 
time, because they are directly proportional to y and y , respectively. For the case without any turbulence, the drop 
distortion motion, according to the breakup model, is created by the difference of forces between the aerodynamics 
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and the sum of the surface tension and viscosity. Obviously, this process is accelerated when turbulence is 
considered. 

A comparison of y and y computed from TAB and T-TAB are presented in the Figure 9. When these forces act 
on the drop, the distortion displacement y would increase with time. Moreover, the imbalance of these forces would 
lead to an increase in the time derivative of this displacement y . When y reaches the value of one, the parent drop 
produces into smaller drops. Figure 9 also reflects a consistency with the results observed in Figure 8. Because of 
an additional turbulence term considered in the T-TAB model, it predicts higher values of y and y at a given time in 

comparison to the results of TAB. Subsequently, the breakup time determined from the T-TAB model is shorter 
than the one from the TAB model. 
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Although no measured data of the secondary breakup are available for the test cases listed in Table 1, they are 
still computed with both the TAB and T-TAB models to assess the variation in their results. Figure 10 shows a 
comparison of the predictions from TAB and T-TAB. The x-coordinate of the plot shows the ratio of aerodynamic 
and turbulence forces. This parameter is selected for plotting since the aerodynamic force is generally dominant in 
the droplet breakup process. The right and left vertical coordinates present the percentage differences in the product 
drop size and breakup time, respectively, predicted from TAB and the present model. The negative scales of these 
two coordinates indicate that the present model predicts smaller values. The results also suggest that the large 
variation of the two predictions appears at the low aerodynamic/turbulence force ratio where the turbulence action 
becomes significant in comparison to the aerodynamic force. However, the variation becomes smaller with an 
increase of this ratio value where the aerodynamic force is dominant. This trend suggests that the large variation 
between the two predictions results from the large turbulence force. When the aerodynamic force is significantly 
higher then it plays a key role in the droplet breakup. For this condition, the results from TAB and T-TAB are not 
much different regardless of how high the value of the turbulence force is. 
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Figure 10. Difference in product drop side and breakup time predicted by TAB and T-TAB 


IV. CONCLUSIONS 

We have presented two new models for modeling turbulence effect in atomizing liquid spray. In the T-blob 
model for primaiy breakup, both the characteristic length and time scales of the surface wave perturbation and the 
ones of the turbulence motion are combined in such a way that their contribution to the breakup mechanism is 
weighted by means of the kinetic energy. It has been observed that the initial turbulence quantities play a key role 
on the jet disintegration. Their values are estimated from the geometry and flow conditions of the injection nozzle. 
For the secondary droplet breakup event, an additional force term, composed of the deformation rate of the surface 
distortion, the dissipation rate of the turbulent kinetic energy, is incorporated into the T-TAB model. Assessments 
of these two proposed atomization models have been performed and several observations and conclusions can be 
drawn as follows: 

1 . In the primary atomization regime, the predictions generally show that the intact core length of the turbulent 
liquid jet is slightly shorter than without the turbulence consideration. In fact, at least in the test cases 
considered in this study, the combined characteristic scales used in the T-blob model result in a shorter jet 
breakup time, leading to a shorter intact core as compared to the blob model. 

2. When the turbulent liquid jet disintegrates and then forms droplets, these drops are generally larger than the 
drops produced by a non-turbulent jet The computational results of the present models reflect this finding, 
particularly in test case H-3, where the product drop size predicted from the T-blob model is as much as one 
order of magnitude larger than the size predicted from the blob model. 

3. In general, the turbulence inside the liquid jet affects the primary breakup of the liquid jet. This turbulence is 
characterized by its fluid properties, gas flow conditions, and initial turbulence quantities, ko and e o . In 

turn, these initial quantities are dependent on the flow conditions and the geometry of the injection nozzle. 

4. According to the results from the T-TAB model, the turbulence inside the parent drop also plays a role in the 
droplet surface distortion. Similar to the primary breakup, the level of the turbulence effect on the 
secondary droplet breakup process depends much on its initial turbulence values. However, this turbulence 
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effect is diminished after a few droplet breakup cycles [21]. This is due to the turbulence reduction from 
one generation of drops to the next, which is purely based on the current assumptions of estimating the 
initial turbulence quantities. 

5. Close examinations of the predicted results for several test cases suggest that the values of the forces 
participating in the droplet deformation vary significantly at the initial time. The aerodynamic and 
turbulence forces are the strongest among them. However, the turbulence force approaches a same order of 
magnitude as the surface tension and viscous forces at the breakup time, while the magnitude of the 
aerodynamic force is still maintained. This observation reconfirms that the turbulence has a considerable 
influence on the secondary breakup, when compared with the surface tension and viscosity effects. 
However, the aerodynamic force also still plays a dominant role in the breakup process. 

6. In contrast to the characteristics of liquid jet disintegration, the secondary droplet breakup mechanism 
produces small drops with a short breakup time when turbulence is considered. The turbulence force term in 
the presoit model is constructed to promote the droplet surface distortion. That formulation leads to predict 
a smaller drop size and shorter breakup time than the results obtained from the existing TAB model. 

This research provides a basic framework to enhance the atomization models and the two new models presented 
in this paper have been incorporated into CFD code for full-field atomization simulation. These results will form the 
basis of a future communication. 
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